%--------------------------------------------------------------------------
% computes the neutral stability curves associated with the 'frozen time'
% eigenvalue problem using a 'weighted' pseudoarclength continuation
% scheme. essentially, the Marangoni number is premultiplied by delta^2 in
% order to make it the same order of magnitude as the wavenumbers. this
% just increases the efficiency of the method.
%--------------------------------------------------------------------------


function [k_c, Ma_c] = weighted_neutral_stab_curves(k0)

% N = 200;
N = 50;
N_search = 30;
ds = 1;

tol = 1e-8;

% preallocate
Ms_c = zeros(N, 1);
k_c = Ms_c;

ev_hat = zeros(N_search, 1);
k_hat = ev_hat;
Ms_hat = ev_hat;
theta = ev_hat;

p = params;


% first and second computation using secant method
for n = 1:2
    
    if n == 2
        p.Ma = 0.90 * Ms_c(1) / p.delta^2;
    end
    
    for its = 1:N_search
        
        if (its == 1)
            if (n == 1)
                k_hat(its) = k0;
            else
                k_hat = k_hat(n-1);
            end
        elseif (its == 2)
            k_hat(its) = k_hat(its-1) + 1e-5;
        else
            k_hat(its) = k_hat(its - 1) - ev_hat(its-1) / dev;
        end
      
        ev_hat(its) = comp_eigs(k_hat(its), p);
        
        if (abs(ev_hat(its)) < tol)
            k_c(n) = k_hat(its);
            Ms_c(n) = p.delta^2 * p.Ma;
            break
            
        else
            
            if (its > 1)
                dev = (ev_hat(its) - ev_hat(its-1)) / (k_hat(its) - k_hat(its-1));
            end
            
        end        
        
    end
end

conv = 1;
n = n + 1;
while (n <= N)
    
    if (conv == 1)
        ds = min(10 * ds, 1);
    else
        ds = ds / 10;
    end
    
    fprintf('n = %d\n', n);
    
    t_vec = [k_c(n-1) - k_c(n-2), Ms_c(n-1) - Ms_c(n-2)];
%     t_vec = t_vec / norm(t_vec);
    n_vec = [t_vec(2), -t_vec(1)] / norm(t_vec);
    
    k_hat(1) = k_c(n-1) + ds * t_vec(1);
    Ms_hat(1) = Ms_c(n-1) + ds * t_vec(2);
    
    p.Ma = Ms_hat(1) / p.delta^2;
    
    theta(1) = 0;
    ev_hat(1) = comp_eigs(k_hat(1), p);
    if abs(ev_hat(1)) < tol
        
        k_c(n) = k_hat(1);
        Ms_c(n) = Ms_hat(1);
        conv = 1;
        n = n + 1;
        
    else
        
        for its = 2:N_search
            
            if (its == 2)
                theta(2) = 1e-5;
            else
                theta(its) = theta(its-1) - ev_hat(its-1) / dev_dtheta;
            end
            
            k_hat(its) = k_hat(1) + theta(its) * n_vec(1);
            Ms_hat(its) = Ms_hat(1) + theta(its) * n_vec(2);
            
            p.Ma = Ms_hat(its) / p.delta^2;
            ev_hat(its) = comp_eigs(k_hat(its), p);
            
            fprintf('%.4e\n', abs(ev_hat(its)));
            
            if (abs(ev_hat(its)) < tol)
                k_c(n) = k_hat(its);
                Ms_c(n) = Ms_hat(its);
                fprintf('-----------------------\n');
                conv = 1;
                n = n + 1;
                break
            end
            
            dev_dtheta = (ev_hat(its) - ev_hat(its-1)) / (theta(its) - theta(its-1));
            
        end
        
        if its == N_search
            fprintf('no convergence!\n');
            conv = 0;
        end
        
    end
    
end

Ma_c = Ms_c / p.delta^2;


small_k = 80 * (1 / p.delta + (1 - 2 * p.beta) ./ k_c.^2) / p.beta / (1 - p.beta);
big_k = 8 * k_c.^3 .* sinh(k_c) ./ (k_c .* cosh(k_c) - sinh(k_c) + 1/3 * k_c.^3 .* exp(-k_c)) / p.delta / p.beta / (1 - p.beta);

kt = logspace(log10(k_c(1)), log10(k_c(end)), 30);
for n = 1:length(kt)
    k = kt(n);
    A = w_coeffs(1, 1, kt(n), p) / p.Ma;
%     an(n) = 0.12e2 * k ^ 4 * sinh(k) / ((-0.2e1 * A(1) * k ^ 3 + (0.3e1 * A(4) + 0.3e1 * A(3)) * k ^ 2 - 0.3e1 * k * A(1) + 0.3e1 * A(3)) * sinh(k) ^ 2 + 0.6e1 * cosh(k) * (A(1) * k ^ 2 + (-A(3) - A(4) / 0.2e1) * k + A(1) / 0.2e1) * sinh(k) + 0.2e1 * cosh(k) ^ 2 * (A(1) * k ^ 2 + (0.3e1 / 0.2e1 * A(4) + 0.3e1 / 0.2e1 * A(3)) * k - 0.3e1 / 0.2e1 * A(1)) * k);
    an(n) = 0.12e2 * k ^ 3 * (p.delta * (1 - 2 * p.beta) * cosh(k) + sinh(k) * k) / ((0.2e1 * A(1) * k ^ 3 + (0.3e1 * A(4) + 0.3e1 * A(3)) * k ^ 2 - 0.3e1 * k * A(1)) * cosh(k) ^ 2 + 0.6e1 * cosh(k) * (A(1) * k ^ 2 + (-A(3) - A(4) / 0.2e1) * k + A(1) / 0.2e1) * sinh(k) - 0.2e1 * sinh(k) ^ 2 * (A(1) * k ^ 3 + (-0.3e1 / 0.2e1 * A(4) - 0.3e1 / 0.2e1 * A(3)) * k ^ 2 + 0.3e1 / 0.2e1 * k * A(1) - 0.3e1 / 0.2e1 * A(3)));
end

figure
loglog(k_c, Ma_c, 'k', k_c, small_k, 'r--', k_c, big_k, 'b-.', kt, an / p.beta / (1 - p.beta) / p.delta, 'ko','linewidth', 1);
% ylim([1e4, 1e6]);
xlabel('$k$', 'interpreter','latex','fontsize',12);
ylabel('$\mathrm{Ma}$', 'interpreter','latex','fontsize', 12);
l = legend('numerical','small $k$ asymptotic','large $k$ asymptotic', 'location','best');
set(l, 'interpreter','latex','fontsize',11);